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A Three Dimensional Gravitational Billiard in a Cone 


Cameron K. Langeij^] and Bruce N. MilleijT] 

Department of Physics and Astronomy, Texas Christian University 

Billiard systems offer a simple setting to study regular and chaotic dynamics. Gravitational 
billiards are generalizations of these classical billiards which are amenable to both analytical 
and experimental investigations. Most previous work on gravitational billiards has been 
concerned with two dimensional boundaries. In particular the case of linear boundaries, 
also known as the wedge billiard, has been widely studied. In this work, we introduce 
a three dimensional version of the wedge; that is, we study the nonlinear dynamics of a 
billiard in a constant gravitational field colliding elastically with a linear cone of half angle 
0. We derive a two-dimensional Poincare map with two parameters, the half angle of the 
cone and £, the ^-component of the billiard’s angular momentum. Although this map is 
sufficient to determine the future motion of the billiard, the three-dimensional nature of the 
physical trajectory means that a periodic orbit of the mapping does not always correspond 
to a periodic trajectory in coordinate space. We demonstrate several integrable cases of the 
parameter values, and analytically compute the system’s fixed point, analyzing the stability 
of this orbit as a function of the parameters as well as its relation to the physical trajectory 
of the billiard. Next, we explore the phase space of the system numerically. We find that for 
small values of £ the conic billiard exhibits behavior characteristic of two-degree-of-freedom 
Hamiltonian systems with a discontinuity, and the dynamics is qualitatively similar to that of 
the wedge billiard, although the correspondence is not exact. As we increase £ the dynamics 
becomes on the whole less chaotic, and the correspondence with the wedge billiard is lost. 
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I. INTRODUCTION 

Billiard systems are a class of simple, analytically tractable models exhibiting the fundamen¬ 
tal aspects of nonlinear dynamics. These systems, first introduced by Birkhoff [1], consist of a 
point mass (the “billiard” or “particle”) in a two-dimensional (convex) region with a piecewise 
smooth boundary, where the motion between collisions is inertial and collisions with the boundary 
are specular and elastic. The natural way to study these systems is using a Poincare surface of 
section taken at encounters with the boundary; this reduces the dynamics to a discrete mapping. 
For different boundaries the motion can be regular, chaotic, or mixed with KAM islands [2]. For 
example, the stadium billiard has been proven by Bunimovich sin to be ergodic while the elliptic 
billiard, first introduced by Berry [2], is integrable. In addition quantum versions of these systems 
have also been studied by Waalkens et al. [5J among others. While these classical billiards are 
optimal for analytical study, more experimentally approachable models accounting for the Earth’s 
gravitational field, called gravitational billiards , have also been widely studied [SHE]. I n [6j, the 
wedge billiard , consisting of a particle falling between two symmetric linear boundaries of angle 
26 , was shown by Lehtihet and Miller to exhibit the full range of possible behavior in Hamiltonian 
systems with two degrees of freedom. Namely, for 6 < 45° the phase space consists of a mixed 
phase space with regular and chaotic regions, for 9 — 45° the system is integrable, and for 6 > 45° 
the system is ergodic. Subsequent work on the wedge billiard includes that of Richter et al. [7], 
where the oscillations in the relative amount of chaotic versus regular parts of the phase space 
(the so-called “breathing chaos”) are discussed in terms of the symmetry lines of the system, as 
well as Wojtkowski [8] who rigorously established the ergodicity of the wedge for 9 > 45°. The 
one-dimensional motion of two particles in a constant gravitational field, a system which is isomor¬ 
phic to the wedge, was studied by Whelan et al. |9j in terms of the stability of fixed point orbits. 
Korsch and Lang m demonstrated the integrability of the motion of a gravitational billiard in 
a parabola, and the corresponding hyperbolic system was investigated by Ferguson et al. in m, 
where the wedge and parabolic behavior was shown to arise in two limiting cases of the parameter 
values. The quantum version of the wedge was discussed by Szeredi and Goodings PHE] in terms 
of Gutzwiller’s periodic orbit theory. More recently, the circular, elliptic, and oval gravitational 
billiards were studied by da Costa et al. m, who showed that the energy of these systems plays 
a key role in separating regular and apparently ergodic behavior. For the case of elastic collisions, 
the possible dynamical behavior in two-dimensional Hamiltonian billiards is well established. 

Classical billiard systems have also been studied in three (and higher) dimensions. The inte¬ 
grable class of ellipsoidal and related billiards were introduced by Richter et al. PS] via action 
integrals, while the semiclassical and quantum versions were studied by Waalkens et al. HZj. How¬ 
ever, less is known about three dimensional gravitational billiards, as these are nonintegrable in 
general. 

One way to study billiards (both classical and quantum) experimentally is by bouncing ultracold 
atoms off of beams of light. These so-called “optical billiards” were introduced by Raizen et al. m 
and provide a testing ground for undriven gravitational billiards. Of course, in any macroscopic 
billiard energy is no longer conserved, so in order to observe nontrivial experimental behavior 
the driven versions of these two-dimensional gravitational billiards must be considered. Recently, 
the driven wedge, parabolic, and hyperbolic gravitational billiards were studied experimentally 
in [IS]. In their experiment, Feldt and Olafsen used a steel ball moving within an aluminum 
wedge, parabolic, or hyperbolic boundary, bounded from above by another aluminum boundary. 
The boundary was driven sinusoidally in the horizontal direction to counteract energy loss due to 
collisions. The results of this experiment showed that the driven parabolic system was regular, 
the wedge chaotic, and the hyperbolic mixed, sharing (as in the static case) characteristics of the 
wedge and parabolic systems at different parameter values. 
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One-dimensional driven gravitational systems have been thoroughly investigated [20H25] . The 
seminal example is the so-called gravitational bouncer , consisting of a particle impacting a period¬ 
ically driven wall in the presence of a constant gravitational field. The gravitational bouncer was 
introduced as a variant of the well-known Fermi-Ulam model H2 E3, and has been studied for 
several types of driving motions including sinusoidal [201 - 122] and piecewise linear [23H25] . To model 
the Feldt experiments, two-dimensional driven gravitational billiards were studied numerically in 
HUES]; in the former theoretical study rotational effects were ignored, while in the latter rotational 
effects were included in the model, making analytical computations e.g., periodic orbits difficult. 
However, in the work of Hartl et al. m, the numerical simulations supported the experimental 
results of HE], except in secondary quantities derived from the data, such as the tangential velocity 
of the particle after collisions. 

One possible issue with experiments conducted on the driven wedge is the two-dimensional 
nature of the idealized system. In the real system, the billiard is not a point particle and thus 
has rotational properties which are affected at each collision. In the experiments of Feldt et al., 
additional boundaries were used to ensure the motion of the billiard was contained in the plane; 
however, any collision with these boundaries would likely play a nontrivial role in the dynamics. 
One way of eliminating this problem would be to get rid of the constraint that the motion be con¬ 
tained in a plane. If, instead of a particle in a wedge, we considered a particle in a cone, then there 
would be no need for additional boundaries. In fact, such a system would be ideal for studying 
the effects of rotation on billiard systems, as the cone could either be driven in the conventional 
sense (i.e., the entire cone oscillating in a fixed direction) or the cone could spin in a sinusoidal 
fashion. As the equations determining the time of the next collision are in a sense unaffected by 
this “rotational” driving, such a system is analytically tractable, while also being experimentally 
realizable. 

In this paper, we take the first step to understanding the conic billiard system. Namely, we 
consider the motion of a uniformly accelerated particle in M 3 , colliding elastically with a linear 
cone of half-angle 9. Along with the energy, the ^-component of the particle’s angular momentum 
is conserved. Using these two integrals of the motion together with the fact that the dynamics 
depends only on the difference in the azimuthal angle <j) between collisions, the conic billiard is 
reduced to a two-dimensional Poincare surface of section with two parameters. However, in order 
to determine the physical trajectory of the billiard uniquely the corresponding change in azimuthal 
angle must be accounted for as well. 

This paper is organized as follows. In section [TTJ we introduce the model and derive the discrete 
mapping characterizing it. In section HI we demonstrate some simple properties of the mapping, 
including integrable limits and the existence of periodic orbits. We analyze the linear stability of 
these periodic orbits in terms of the parameter values. In section IV we present some numerical 
results, and in section[V]we discuss the results obtained, as well as possible extensions of our model. 


II. THE MODEL AND THE MAPPING 


As the motion of the particle is unaffected by its mass, we are free to set the mass of the particle 
equal to unity, without loss of generality. We orient our Cartesian coordinate system so that the 
axis of the cone is the z- axis, and 9 is the usual spherical polar angle [29] . The gravitational field is 
g = —ge z , where e z is a unit vector in the z-direction. We shall find that the Poincare map is most 
easily obtained in spherical polar coordinates (r, 0, </>). In this coordinate system, the Hamiltonian 
of the billiard between collisions is 





+ gr cos 9. 


(i) 
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This governs the motion of the particle between collisions, which we assume to be elastic. There 
are two independent integrals of the motion preserved by both the parabolic motion and collisions 
between the particle and the cone: 


H = E, and p^ = l z , (2) 

which we identify as the energy and z-component of angular momentum, respectively. The presence 
of two conserved quantities reduces the dimension of the phase space of the system from six to 
four. Taking the natural Poincare surface-of-sect ion at the moment of collision (here p — \[x 2 + y 2 
is the cylindrical polar coordinate), 

z = p cot 6 = r cos 6 (3) 


further reduces the dynamics to a three-dimensional map. In fact, as we shall show below, the 
azimuthal angle 0 plays no role in the dynamics, and therefore the arbitrary initial condition 
0 = 00 is sufficient to describe all possible behavior of the system. Thus, a two-dimensional map 
characterizes the dynamics. By a suitable transformation of the coordinates and time, the constants 
E and g can be re-scaled arbitrarily; for convenience, we choose units such that E = g = so 
that the three constants E,g,£ z are consolidated into the single (dimensionless) parameter (where 
the factor of ^ is chosen for aesthetic purposes) 


n! — 

V2EV 2 ’ 


( 4 ) 


With this choice of units, energy and z-component of angular momentum conservation are expressed 
by (where we have replaced the momenta p r ,po with the corresponding velocities v r , vq) 


E 2 

1 =v1+Vq+ , . 9 ^ + r cos 0 , E = rvs sin 6. 
r z sin 0 


( 5 ) 


The cone angle 6 is restricted to 0° ^ 6 ^ 90°. It is not difficult to show using energy conservation 
that E satisfies \E\ ^ tan 6. In fact, a more careful analysis shows that this bound can be lowered 
to \E\ ^ calling this upper bound £ max , we obtain a parameter which ranges from zero to 

one by defining £ = -/—■ This defines the system and its parameters. 


We have shown that a two dimensional map is sufficient to characterize the dynamics of the 
conic billiard. However, we not yet shown which variables are most suitable. It is clear that 
this choice cannot be made haphazardly, as most pairs of conjugate variables are insufficient to 
fully determine the particle’s subsequent trajectory. As we shall see below, it turns out that the 
radial coordinate and the radial component of the velocity, (r, v r ) provide enough information to 
determine the subsequent motion unambiguously (for a given choice of initial 0). Thus, our aim is 
to compute the Poincare map 


V : 



( 6 ) 


where the suffix n denotes the value of a quantity just after the nth collision. In order to compute 
this map, we require two ingredients: (1) the time interval between collisions, and (2) the law 
relating the pre-collision and post-collision components of the particle’s velocity. To compute (1), 
we need only basic kinematics; in Cartesian coordinates, the time of the (n + l)st collision is defined 
by the implicit relations 


Vn x (fn+l tn) T X n — 3?n+l 


Vriyifn+l tn) T Un — Un+l 5 


( 7 ) 
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- \ttn +1 ~ tnf + V nz {tn +1 “ t«) + \/®n + Cot 0 = ^/*n+l + ^n+1 COt0 ' ( 8 ) 

where the factor | is due to the fact that we set g — and we have used the relation 2 : = pcot 6 
defining the surface of the cone. These equations are equivalent to the cubic equation (where 
T~n+ 1 = Ui+1 — tn) 


T n+1 ~ 8 V n ,T^i + 16 


J n z 1 n+1 


<-«+»!)cot 2 0- v'^ + y. 2 . 


cot 0 


^n+l 


+ 32 


V Hz + Vn cot 0 ~ ( x nV nx + VnVriy ) COt 2 0=0, (9) 


which can be expressed as 


x n+l + 8 ( v n e Sm0 - V Ur COS0) 

lT n+l 

+ 16 jn 2 e - 

( v n e + v l *) cot 0 + 


r n cos 6 \ 

~^2~ J Tn+1 

— 32 r n v ne cot 6 = 0 (10) 


in spherical coordinates. The smallest positive root of this cubic equation is the time of the (n+l)st 
collision. The roots of this cubic equation are long and complicated, and do not yield any physical 
insight since a priori there is no easy way to unambiguously determine the smallest positive of the 
three (real) roots. Thus, they are not included here. 

To analyze the effects of a collision, spherical coordinates are natural for this system. This is 
because the equations relating pre-collision and post-collision components of the billiard’s velocity 
become especially simple in these coordinates, since e# is orthogonal to the surface of the cone, and 
e r , are parallel to it (where e r , e#, are unit vectors in the direction of r, 6 and </>, respectively) 
at the moment of collision. If v~ +1 . denotes the ith component of the particle’s velocity just before 
the (n + l)st collision, and u^ +1 . the corresponding component just after the (n + l)st collision, 

then the vector equation (here n = e# is a unit normal to the surface of the cone, and is the 

velocity of the particle just after/before the (n + l)st collision) 

V+ +1 = v“ +1 - 2(v“ +1 -n)n (11) 

is equivalent to the three component equations 


V n+l r V n+ l r > V n+ 1^ V n+ 1^’ V n+lo V ti+Iq' (^) 

Now we demonstrate that knowledge of r and v r at the nth collision is sufficient to describe the 
subsequent motion of the particle. Suppose these values, (r n ,u nr ) at the nth collision are known. 
Then we immediately know the value of 6 n = 0, since at each collision this is just the half angle of 
the cone; additionally, we know the ^-component of the particle’s velocity, from the conservation 
of the z-component of angular momentum: 


v 


_ £ 

U(t) r n sin 6 ’ 


(13) 


Using this in the energy conservation expression permits us to solve for in terms of r n and v Ur : 


v 


2 

no 


1 — r n cos 6 — v 


2 

n r 


i 2 

sin 2 0 


(14) 
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Although it seems we may only know v ne up to a sign, in fact we are safe to choose the negative 
root, since v$ is the component of the velocity orthogonal to the surface of the cone in the direction 
of e#, which, due to the convention e# = x e r , cannot be positive after a collision. Moreover, 
as we shall see below, the azimuthal angle <f n plays no role in the dynamics; thus, it can be fixed 
arbitrarily at t — 0 and subsequently eliminated from the mapping, as only the difference (j) n +i ~ <fn 
appears in the mapping equations (and this quantity can be computed purely in terms of r n ,v Ur 
and the time of the next collision). Thus, to complete the mapping we need only compute the 

time interval from (10), the new radial distance r n+ i = yj x 2 +1 + p 2 +1 csc $ an d the radial velocity 
component v n +i r from (see Appendix) 


v n +i r = v nr (sin 2 0 cos p n +i + cos 2 6) + v nQ (cos tp n + l - 1) sin 6 cos 6 


+ v n4) sin (p n+ i sin 6 - -r n+ i cos 6 , (15) 


where ip n +i = (j) n + 1 — </>n is the difference in azimuthal angle between collisions, given by 


<Pn+i — arctan 


v n^ T n J r 1 

r n sin 9 + (v Ur sin 9 + v ne cos 9)r n+1 J ’ 


(16) 


with care being taken to choose the correct quadrant in the xy- plane. In terms of these quantities 
the mapping for r n+ \ is 


1 


2 

(' V nr + V ne Cot 0) + ~2 


T% +1 + 2 r n {v Ur + v ne cot 9) T n+ 1 + rl, 


(17) 


and v n +i r is given by (15). After evaluating sin(/? n+ i and cos(/? n+ i using (16) and defining the 
convenient reduced variables p = rsin0, u r = v r sin 9, uq = vqcos9, the mapping can be rewritten 
as (see Appendix) 


2 

Pn +1 


( U nr + U, 


n e j 


e 

H-2 

Pn J 


r 2 +1 + 2 p n (u nr + u ne ) r n +i + p, 


2 

Vi 5 


(18) 


l r 


sin 2 9 
Pn+l 


{u. 


2 (f~\ 1 

+ u n e ) H 2 Tri +1 + Pn ( u n r + u n e ) r + u n r C °S 2 0 ~ 

Pn _ J 


m r 




Of course, our ‘reduced variable’ p is just the cylindrical polar coordinate; however, u r is not the 
corresponding velocity component v p . In fact, v p = u r + uq. Although this suggests that the 
conjugate variables (p, v p ) should be used for the mapping, these two variables are insufficient to 
determine the subsequent trajectory of the particle. This can be seen by writing energy conservation 
in terms of cylindrical polar coordinates: 


1 


£ 2 

= v 2 + v 2 H— 2 T p cot 9. 


Thus if we know p and v pi the square of v z is given by 


2 n 2 
v z^ 1 ~ v p 


t 2 

— pcot 9. 
P 2 


( 20 ) 


( 21 ) 
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However, without some additional piece of information regarding the sign of v z , a priori we have 
no way of determining which square root to take. Thus (r, v r ) or equivalently (p, u r ) are the 
coordinates we use in this work. Although in principle we can eliminate vq or u$ from the mapping 
via energy conservation, the equations become considerably more complicated and yield little 
additional insight. Moreover, we find that computing periodic orbits of the system is easier when 
using this ‘implicit’ form of the mapping. 

The reason we choose (r, v r ) as our primary mapping variables is simple: they are area-preserving 
i.e., the determinant of the Jacobian matrix J = - g e q Ua j un ity. Another sensible 

0{r n ,Vn r ) 

choice would be (sgn(v r )v\\,v±), where and v± are the tangential and normal components of the 
velocity with respect to the cone boundary, respectively. Although these variables are sufficient to 
fully determine the dynamics, they are not area-preserving. In the next section, we shall utilize 
several times the area-preserving nature of (r,v r ) in the computation of stability eigenvalues. In 
deriving specific fixed points and periodic orbits, the reduced mapping (18)-(19) is more suitable. 
Since the variables are not area-preserving, however, they are not ideal for surface-of-section plots. 


III. ANALYTICAL RESULTS 


In this section we examine some general properties of the mapping. In section IIIA 


we examine 


limiting cases of the system’s parameter values, demonstrating that the conic billiard becomes 


integrable in several limits. In section [Til B| we compute the system’s fixed points and analyze their 
linear stability as a function of parameter values. 


A. Simple properties of the mapping 


We note that if £ = 0 the conic system evidently reduces to the two-dimensional wedge billiard, 
as the Hamiltonian becomes 




wedge 


p 2 r + ) + r cos 9. 


( 22 ) 


Hence, following [b] we see that for t — 0 the system becomes integrable for three cases (i) 0 —>> 0°, 
(ii) 9 — 45°, and (iii) 9 90°. In the first case, the potential becomes purely radial in the limit 

i.e., 


H 


Pr 


+ 


P 2 e 


+ r , 


(23) 


implying the conservation of p$ = r 2 0, which is also preserved by the collision map. For 9 = 45°, co¬ 
ordinates parallel and orthogonal to the wedge surface can be defined so that the motion becomes 
separable. Finally, for 6 = 90° the motion becomes simply that of a projectile bounded from below 
by a horizontal floor. In this case, the motion is unbounded. 

For £ / 0, we lose integrability at 9 = 45°. However, we retain the integrable limit 9 —» 90°, 
as the above argument is unaffected by the introduction of a nonzero z-component of angular 
momentum. Additionally, the limit 0—^0° remains integrable as long as we simultaneously take 
t ^max- Examining the Hamiltonian 





+ 


e 


r 2 sin 2 6 


+ rcosO, 


(24) 









we see that if we naively take the limit 0 —>► 0 the term proportional to sin 2 0 diverges. However, 
if we instead take 0 —> 0 with £ —»> ^ max ^ tan 0 , we see that 

H 2 if r + ) + ^ + r ’ ( 25 ) 

where a = Hence the potential again becomes indistinguishable from a central potential in this 
limit, implying that the quantity p@ is conserved and the motion integrable. 


B. Fixed points 


Imposing the conditions p n +i = p n = p, ^n+i r = u nr = on the mapping, we arrive at (here 
r n+ i — r — const.) 


0 = 


(u r + Uq) 2 H—« 


T + 2 p(i/ r + TZfl), 


(26) 


sin 2 0 


u r = 


*2 1 


(?i r + ^6>) H-o 

p J 


r + p(u r + uq) > + u r cos 2 9 — uq sin 2 9 — — sin 29. (27) 


The first equation is equivalent to 


T = 


2 p(uj> 'Uq') 

(u r + ^ 6>) 2 + £ 2 1P 2 ’ 


(28) 


Evidently this holds for all fixed points of the map. Using this in the mapping equations for u n +± r 
and u n +i e gives the simultaneous equations 


r = —4 {u r + uq) tan 9 


and 


r = —4 u r cot 9 — \uq tan 9. 


(29) 


Thus either u r — 0 or 9 = ?. In the case u r — 0, setting the two expressions for r equal gives 


£ 2 

p — 2 tan 9 ( Uq H—^ 


If f = 0, then it is easily verified that 


P = 


sin 20 
2 + cos 29 1 


u r — 0 


(30) 


(31) 


is a fixed point; this agrees with the fixed point of the wedge billiard [ 6 ]. In the general case that 
£ 7 ^ 0 we can apply energy conservation to arrive at 


Uq 


(p 2 — 3£ 2 ) cos 2 9 
p 2 (2 + cos 20 ) 


(32) 


which must hold simultaneously with ([30]). Eliminating uq in ([30]) and (32) yields the cubic equation 

[(2 + cos 29) cot 9} p 3 — (2 cos 2 0) p 2 — 2£ 2 sin 2 0 = 0. (33) 

Clearly, if we set £ — 0 in the above equations we recover the limiting case found above. Although 
the equation determining p is cubic, it turns out that for all allowable parameter values only one 
of the three roots is real. Hence there is a single fixed point 


p 


u r = 0 


where p* solves (33). 


(34) 
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The stability of this orbit can be analyzed in terms of the eigenvalues of the Jacobian matrix. As 
is shown in e.g., [2], the stability condition that the eigenvalues of the Jacobian are less than one 
in modulus is equivalent to, for an area-preserving map, the condition 


|TrJ| < 2. 

This can be alternatively stated in terms of Green’s residue [3D], defined as 

R= 2 -^. 


(35) 


(36) 


In terms of i?, a fixed point is stable for 0 < R < 1, and unstable for R < 0 and R > 1. 
Stable fixed points are also called elliptic, and unstable fixed points are called hyperbolic. Of 
course, since the (p,u r ) map is not area-preserving, we must use the Jacobian of the (r,v r ) map 
in stability calculations. In terms of these variables the fixed point is 


r = r*, v r = 0, 

where r* is given by 

[(2 + cos 26) cos 6] r* — (2 cos 2 6) r 2 — 2 1 2 — 0. 
The real root of this cubic can be written 

4 cos 2 6 2 cos 9 


r* = k(£, 6) + 


+ 


9(2 + cos 2 9) 2 k(£, 9) 3(2 + cos 29) ’ 


(37) 


(38) 


(39) 


where k(£,9) is 


k(e,e) = 


8 cos 3 9 


+ 


£ 2 sec 9 


+ 


7 sec# 


_27(2 + cos2 9) 3 2 +cos 2 9 3+3(2 + cos 29) 2 


- vT6a3s4?T27F(2Tcos20)2 


1/3 


(40) 

It is easy to see that this agrees with the case £ = 0 in equation (31). Since the mapping depends 
explicitly on the time between collisions r n+ i, in computing the necessary derivatives we must take 
care to use the chain rule i.e., 


dr 


71+1 


dn 


77+1 


d r„ 


dr„ 


<9r n+ i dr n+ i du n+ i r _ <9u n+ i r 

r n+ i=const. dr n+ i dr n ’ dv Hr dv nr 


+ 


+ (41) 


r„+i=const. dr n - |_i dv 


n r 


where the derivatives of r n+ i are computed using implicit differentiation on the cubic equation for 
the time interval i.e., 


dTn+l 

dr n 


df/drn 

df/dr n+ i 


(42) 


where /(r n +i, r n , v Ur ) — 0 is defined by the left hand side of ©■ The computations are straight¬ 
forward but tedious due to the fact that v ne depends on r n and v Ur in the energy expression (14). 
The explicit expression for Tr J, or equivalently Green’s residue R is long and complicated, and 
will not be included here. Instead, in Fig. [l](a) we display the stability regions in (£,9)- space. 

We note that although (r*,0) is a fixed point of the mapping , the physical trajectory of the 
billiard does not in general form a closed orbit. Instead, fixing r = r* at each collision (which is 
equivalent to fixing the height of the particle at each collision) constrains the collision points to 
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FIG. 1. (a) Stability regions of the fixed point in the parameter (£, 0)- space. Black regions correspond to an 
elliptic (stable) fixed point. Note that for £ — 0 the stability region is precisely 0 ^ 9 < f. (b) ^-projection 
of the real-time trajectory of the billiard for 9 = ^ and £ = 0.1 for the first 20 collisions of the fixed point. 


a circle on the conic boundary surface. It is clear that the trajectory of the billiard will form a 
closed path if and only if (note that <p n +i = (p = const, for a fixed point) 



where p, 


(43) 


Physically, this orbit corresponds to a trajectory which repeats after q collisions, during which the 
billiard’s ^-coordinate cycles through 2np radians. When p is not a rational multiple of 2 tt, the 
physical trajectory is quasiperiodic. Since the fixed point orbit is confined to constant z and v z , the 
physical trajectory corresponding to this orbit is equivalent (when viewed from the positive z- axis) 
to that of the circular billiard i.e., a billiard moving uniformly within a circle. In Fig.[l](b) we plot 
the projection of the physical trajectory in the xp-plane for 6 — 30° or | radians and £ = 0.1. 

For more general periodic orbits, we note that for £ = 0 there are, for example, well-known [6] 
period-m orbits of the form 


p = tan 6 


1 - 


m 2 cos 2 6 + sin 2 6 
1 — (m + l) 2 cos 2 6 cos 2 6 


m sin 9 cos 9 


^/l — (m + l) 2 cos 2 9 cos 29 


(44) 


Such solutions clearly hold for the conic billiard when £ = 0; however, generalizing these solutions 
to £ ^ 0 is difficult to do analytically due to the implicit nature of the mapping. However, the 
existence of higher order periodic orbits can be confirmed numerically. 


IV. NUMERICAL RESULTS 


In this section, we iterate the Poincare map numerically to examine the qualitative changes in the 
phase space as the parameters of the system are varied for two values of £ which are representative 
of the dynamical behavior. In section [IV A we consider a relatively small value of £ — 0.1, which 
corresponds to small values of the particle’s ^-component of angular momentum. Intuitively, we 
expect the behavior of the system in this case to be qualitatively similar to that of the wedge billiard 
pm. In section [IV B[ we increase the particle’s angular momentum to £ — 0.5 and examine how 
the wedge-like behavior is destroyed and, in the limit £ 1, approaches integrability. 
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FIG. 2. SOS at £ = 0.1 for (a) 9 = 15°; (b) 9 = 18.5°; (c) 9 = 21°; (d) 9 = 24.5°; (e) 9 = 27°; (f) 9 = 30.5°. 


A. Wedge-like behavior: small 


For £ — 0 it is clear that the conic billiard reduces to the wedge billiard. For relatively small 
values of 0 < £ < 0.2 we find that some structures found in the wedge system survive; however, 
the correspondence is not exact. Specifically, the conic billiard does not in general become ergodic 
for angles above an exact value (45°, for the wedge). Instead, for small £ values the conic system 
becomes ergodic in only a certain range of angles, with stable periodic orbits reappearing for large 
angles. 

In Figs. 2]|4 we display the Poincare surface-of-section ( r n ,v nr ) for a variety of cone angles at 
£ = 0.1. We note that the general structure of the phase space is quite typical of two-dimensional 
area-preserving maps. As 9 is varied, we see a number of interesting bifurcation processes as 
periodic orbits lose and gain stability. The fixed point found in section |IIIB 
surrounding KAM islands, is seen in Fig. [^ a 


together with its 
5, 9 < 30° to constitute 


-(f), which corresponds to 15' 
the majority of the phase space along with chaotic regions. Additionally, for this range of angles 
there are a variety of period-2,3 and higher orbits with their own surrounding island structures. 
Note that in these figures the elliptic nature of the fixed point is clear, which agrees with the 
stability calculations of section IIIB In Fig. [3]^a) at 9 — 34° the fixed point becomes hyperbolic, 
and aside from a single period-3 orbit the majority of the phase space becomes chaotic. As 9 is 
increased further to 9 ~ 45° the stability regions grow, leading to bifurcations and higher periodic 
orbits, as seen in Fig. |3^d)-(e). Finally, in Fig.^f) ah periodic orbits disappear except for a single 
period-2 orbit, which remains until 9 ~ 54°, where the system appears to be ergodic. This behavior 
is analogous to that of the wedge. 

The apparent ergodicity remains for 54° < 9 < 73°, where a stable period-2 orbit re-appears, as 
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shown in Fig. 
point of Fig. [4 


4jb). As 9 is increased to 9 ~ 74.5°, the stability islands collapse to the single fixed 
. As seen in Fig.[4](d)-(f), higher periodic orbits gradually appear and the stability 


region of the fixed point grows as the angle of the cone is increased further. Finally, for large cone 
angles the system approaches integrability, with the chaotic part of the phase space vanishing in 
the limit 6 90°. 


B. Large £ values 


As £ is increased, the relative amount of chaos in the phase space becomes smaller. As we 
increase the allowable we lose the “wedge-like” behavior seen in Figs. [2][4] and the dynamics 
become more regular. In Fig. [ 5 ] and Fig. [6] we plot the Poincare surface-of-section ( r n ,v nr ) for 
£ = 0.5 and a variety of 9 values. The value £ = 0.5 is representative of the behavior of the system 
in the range 0.1 < £ < 0.7. As £ is increased beyond this range, the behavior approaches the 
integrable limit £ —» 1, with the amount of chaos becoming negligible (and disappearing completely 
in the limit). Additionally, we find no instances of a simply connected region of chaos. 

For 0 ^ 9 < 10° the system’s behavior is nearly integrable, as may be seen in Fig. [Ja). The 
fixed point found in section IIIB is once again seen to play a significant role in the phase space, 
with the chaotic region remaining quite small until 9 ~ 42°. As 9 is increased further, a period-4 
orbit shown in Fig. [5](d) collides with a stable island surrounding the fixed point, leading to the 
chaotic band seen in Fig. [5^e) at 9 — 44°. This chaotic band remains and reaches its maximum 
size for 9 ~ 63.5°, as seen in Fig. [6](b) . For larger cone angles the chaotic band gradually shrinks, 
vanishing in the integrable limit 9 90°. 
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(g) (h) 


FIG. 4. SOS at £ = 0.1 for (a) (9 = 54°; (b) 0 = 73°; (c) 0 = 74.5°; (d) 6 = 77°; (e) 0 = 80.5°; (f) 0 = 84°; 
(g) (9 = 87.5°; (h) (9 = 89.5°. 


The transition between the small and large Gbehavior is found to occur for the most part in the 
range 0.15 < £ < 0.25. In general the result of increasing £ is to reduce the amount of chaos in the 
phase space. In Fig. [7] we show this effect by plotting the Poincare surface-of-section for different 
£ values at a fixed 6 = 15°. We see that the chaotic region present at £ = 0.1 has almost entirely 
disappeared by £ — 0.25, and has been replaced by stable island structures throughout the phase 
space. 
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r n 


X n 





(a) 


(b) 


( c ) 



FIG. 5. SOS at £ = 0.5 for (a) (9 = 10°; (b) 0 = 25.5°; (c) 0 = 34°; (d) 0 = 42°; (e) (9 = 44°; (f) 0 = 50°. 


V. CONCLUSIONS & OUTLOOK 


In this work we introduced a three-dimensional gravitational billiard consisting of a particle 
falling in a linear cone. We reduced this system to a two-dimensional area-preserving discrete 
map with two parameters. However, as demonstrated in section IIIB the (r,v r ) mapping must 
be supplemented with the corresponding change in </> in order to fully understand the physical 
trajectory of the billiard. Thus, in this sense, the two-dimensional mapping does not always tell 
the whole story. We found several integrable limits of the system, computed the fixed point of 
the map and examined its linear stability as a function of parameter values. Due to the three- 
dimensional nature of the physical trajectory, this fixed point of the map does not in general 
describe a closed orbit in coordinate space. Next we investigated the phase space of the conic 
billiard in terms of the two parameters, 9 and £. For small £ values we found that the system’s 
phase space was qualitatively similar to the wedge, with a mixed phase space for 0° < 0 < #*, which 
gives way to apparent ergodicity for 9 > 9*. In contrast to the wedge, however, this ergodicity 
does not hold for 9* < 9 < 90°; instead, stable periodic orbits re-appear, eventually leading to 
the integrable limit 9 90°. As we increased £ we found the departure from the behavior of the 

wedge billiard to be dramatic. The relative amount of chaos in the phase space is reduced, with 
no simply connected region of chaos appearing for £ > 0.15. The fixed point of the system plays 
a more significant role for large values of £, being stable for all cone angles with £ > 0.4. This 
fixed point, together with surrounding KAM islands and additional periodic orbits, constitutes the 
majority of the phase space, with only a small chaotic band appearing for certain angles. 

In future work on the conic billiard we would like to include rotational effects, as well as examine 
driven versions of this system. One possible way to introduce time-dependence to the conic billiard 
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FIG. 6. SOS at£ = 0.5 for (a) 6 = 60°; (b) 6 = 63.5°; (c) 0 = 67°; (d) 0 = 70.5°; (e) 0 = 77°; (f) 0 = 81.5°. 


would be to ‘spin’ the cone sinusoidally. This would have the advantage that only the collision 
equations would need modification; computing the time of the next collision, which is in general 
a nontrivial numerical task, would be no more difficult than in the undriven conic system studied 
here. The conic billiard therefore offers an opportunity to examine the effects of rotation and 
time dependence on gravitational billiard systems. Another direction for future work is to consider 
different boundary surfaces. In two dimensions the parabolic billiard is the only known example of 
an integrable gravitational billiard m-, in fact, it is possible to show that the corresponding three 
dimensional system, consisting of a gravitational billiard in a paraboloid, is also integrable. We 
plan to investigate this paraboloidal billiard in more detail in subsequent work. 

The conic billiard is distinct from the wedge primarily because of the dependence of the dynamics 
on the ^-component of the particle’s angular momentum. This is in contrast to most billiard 
systems, where the dynamics is governed only by the shape of the boundary i.e., in this case the half¬ 
angle 9 of the cone. The conic billiard is an example of a system where the initial condition of the 
particle (the (j) component of the velocity) is folded into a parameter of the mapping characterizing 
the dynamics. Just as the properties of the wedge billiard have been experimentally confirmed 
using optical billiards m, the conic system studied here could likely be tested using ultracold 
atoms bouncing off laser beams. 


Appendix: Derivation of the velocity map 

Writing the collision equations 

V n+l r — V n+ l r ’ V n+l ( j ) ~ V n+ 1^’ V n+l e ~ ~ V n+ 1 0 ’ 


(A.l) 
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FIG. 7. SOS at 0 = 15° for (a) £ = 0.1; (b) £ = 0.13; (c) £ = 0.16; (d) £ = 0.19; (e) £ = 0.22; (f) £ = 0.25. 


in terms of the Cartesian velocity components results in a linear system of equations for the post¬ 
collision velocities ^n-\-i x ^ v n-\-i y an< ^ v n+i z t erms of the corresponding pre-collision velocities 
v n+i x ’ v n+i y an< ^ v n+i z - The solution of this system can be written 

V n+I x = V n+I x (! - 2 cos2 </>n+i cos 2 0) - v~ +ly sin 2<p n+1 cos 2 9 + v~ +lz cos <f> n+1 sin 29, 

V n+ly = — ^n+lj; sin 2< ?Wl COS 2 # + Vf 1 (1 - 2 sin 2 </>n+l cos2 9 ) + V n+1 2 sin 0n+l sin 20, (A.2) 

u n+l, = u n+l a C0S 0n+l sin 29 + v n+ l y sin 0n+l Sm 20 + V~ + ^ COS 20. 

Since between collisions the billiard follows a simple parabolic trajectory, we have (recall that we 
set g = \ by a scaling transformation) 


J n+ l a 


— 


v n+l y — ^3/ > 


Ui+1 




+ 'Gi 


(A.3) 


Using this in (A.2) and writing the Cartesian components of the velocity in terms of the spherical 
components, the velocity map becomes (where we have abbreviated v^ +1 . = v n +i i , and defined 
(fn+l — 0n+1 — fin) 


v n+ i r = v nr (sin 2 0 cos (f n+ i + cos 2 0) + v ne (cos <p n+1 - 1) sin 0 cos 0 + v n<j> sin ip n+ i sin 0 

-h n+ icos0, (A.4) 


V n +I e = v nr (1 - cos <y?n+i) sin 0 cos 0 - (cos 2 0 cos <^ n+ i + sin 2 0) — v n , sin <p n+ i cos 0 


— 2 Tn+1 sin ^’ ( A - 5 ) 
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v n +4 = ~v nr sin p n+ 1 sin 0 - sin + n+ i cos 9 + v n<j> cos <p n+ i. 

The difference in azimuthal angle <f n +i is easily found to be given by 

Uin / T 77 -L 1 

<Pn +1 = arctan - . -^- 

|_r n sin 6 > + ( v nr sin 9 + v ng cos 9 )r n+ 1 _ 

and after some algebraic manipulation this gives 

Vn^n+l Vn^Tn+l 


sm (fin+l = 


yj[r n sin 9 + (u nr sin (9 + cos(9)r n+ i ] 2 + rn +1 sm 9 


and similarly 


COS (Pn+l — 

In terms of the reduced variables 


r n + ( T OOt 

Gi+1 


p = r sin 0 , = v r sin 0 , uq = cos 0 

the (p, ix r ) mapping becomes 


Pn+1 


2 

(^n r + U ne ) H-y 

Pn. 


T n +1 + ^Pn (+n r + U ng ) + 


(A. 6 ) 


(A.7) 


(A. 8 ) 


(A.9) 


(A. 10 ) 


^7l+l r — 


sin 2 6 
Pn +1 


\2 ^ 


(^n r +^n 0 ) H-y 

Pn _ 


T n+ 1 + p n ( u nr + U ne )\ + u nr cos z 9 - u ng sin " 1 9 


~pTn+i sin 29. (A. 11 ) 


The cubic equation for r n +1 in terms of these reduced variables is 
Tn+i + 8 {u ng tan 9 - u Ur cot 9 )t z +1 


+ 16 


(sec 2 6 — esc 2 9 )u 2 — 


£ 2 cot 2 6 


2 o\ 2 ^ w o 2 n Pn c 0^0 

-2iq, iq^csc 0 - 


ne 2 

rn 


*n r ^uq 


7"n+1 

,2 


— 32 p n u ne esc 2 9 — 0 . (A. 12 ) 
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